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Abstract 

We present a numerical study of a trapped binary Bose-condensed gas by 

solving the corresponding Hartree-Fock equations. The density profile of the 

binary Bose gas is solved with a harmonic trapping potential as a function 

of temperature in two and three dimensions. We find a symmetry breaking 

in the two dimensional case where the two condensates separate. We also 

present a phase diagram in the three dimensional case of the different regions 

where the binary condensate becomes a single condensate and eventually an 

ordinary gas as function of temperature and the interaction strength between 

the atoms. 
03.75.Fi,05.30.Jp 



1 



Typeset using REVT^ 



I. INTRODUCTION 



The recent realization of Bose-Einstein condensation in dilute gases of trapped alkali 
atoms has provided an opportunity to investigate macroscopic quantum effects in novel 
systems The second generation of experiments on condensates in Rb and Na vapors 

give detailed information about the condensate long range order which distinguishes 
them from an ordinary gas. These atoms are known to have a positive scattering length , 
i.e the atoms experience a repulsive force between them. Recent experiments with Li |^ 
show that, for a limited range of particle densities even particles with negative scattering 
length can form a condensate. In a recent experiment one has been able to trap atoms of 
two different spin states and cool them below the condensate transition point [^f. These 
experiments give us knowledge about the condensation in a whole knew dynamical regime. 

Much work on the density distribution of the condensate has been carried out for a single 
condensate and lately also for the two-component condensate |llO|-p!^. These calculations 
have mainly been valid at the temperatures T = 0. In our earlier paper we have solved the 



Hartree-Fock equations for a single condensate [0]. In this paper we derive and solve nu- 
merically the Hartree-Fock equations for a two-component gas trapped in harmonic external 
potentials at a finite temperature. We solve for the density in two and three dimensions for 
the case with simple harmonic potentials. The Hartree-Fock equations give a phase transi- 
tion at the onset of condensation but they do not take into account any critical fluctuations. 
It is however possible to obtain a qualitative picture of the features of the Bose-Einstein 
phase transition using the HF approach. 

In a real experiment, the different potentials experienced by the two components make 
them respond differently to gravity. This separates the condensed clouds in space, but we 
assume this effect to be compensated by technical means. In this way we can investigate 
the intrinsic effects determined by the physically more essential interactions. The two com- 
ponents of the condensate repell each other, and if they form exactly on top of each other, 
the more weakly interacting one is pushed to expand away from the center. This is found to 
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lead to a condensate component that has a minimum at its center. This situation is still not 
optimal, because energy can be lowered by separating the condensates. In a symmetric trap 
this implies a breaking of the symmetry; the two peaks are shifted in a direction not fixed by 
the external potentials. We investigate this symmetry breaking in a two-dimensional model 
trap, because solving for the asymmetric condensate proves to be numerically demanding. 
The calculations are found to verify the intuitive picture described. 

If the two components did not interact, the condensates would form and disappear in- 
dependently, thus showing two different transition temperatures. With condensate interac- 
tions, the presence of one condensate will affect the formation of the other one. This will 
shift the transition temperatures depending on the strength of the interaction between the 
two components. In particular, we expect a repulsive interaction to work against the simul- 
taneous formation of both condensates. Thus we expect to see a modification of the phase 
diagram in the plane of temperature and interaction between the components. In order to 
keep the numerical effort manageable, we investigate this effect in a spherically symmet- 
ric three-dimensional trap. We find that the interaction shifts the boundary between the 
two- and one-condensate regions. We even find a case where there is no formation of a 
single-condensate region. 

The organization of the paper is as follows: Section II reviews the HF equations that 
we solve numerically. In section III we briefly go through the numerical methods that have 
been used. The results of the calculations are presented in Sec. IV where various situations 
are compared. Finally Section V comments on the calculations and their results. 

II. THE HARTREE-FOCK EQUATIONS FOR THE TWO-COMPONENT 

CONDENSATE 

The Hamiltonian for an interacting Bose gas of two different kinds of atoms can be 
written in the form 

H^H^ + H^ + Hs, (1) 
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where 
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H2 = J dr4{r)[-^V^ + U2{r)]Mr) + 1/ dr J drVl(r')4(r)T4(r - r')V'2(r)V'2(r') (2) 
H3 = ^JdrJ rfrVl(r')4(r)V^,,i(r - r')^i(r)^2(r') 

and "01, ^/^i, ip2 and il^l are the Boson field operators which obey the commutation rule 

[^,(r),^](r')] = (5,,(5(r-r'). (3) 

Ui{r) and Ui{r) are the two different external traps that confine the atoms. We are here 
going to use the short range approximation both for the interaction between the particles of 
the same kind and between the two different kinds of particles. This gives us the interaction 
potentials 

Vi{r - r') ^ vi5(r - r') (4) 
V2{r - r') = V26{r - r') (5) 
VUr-r')^w5{r-r'), (6) 

where 

4:Tih^ai A'Kh'^a2 47r^^ai2 

Vi = V2 = w = [7) 

mi m2 ^ymlm2 

with the s-wave scattering lengths oi, 02 and 012. Our goal is to calculate the temperature 
dependence of the two gases. We therefore introduce the thermodynamic free energy fl{T, /i) 
which is defined as 

g-/3n^y^[g-/3(H-M.ivW)]^ (8) 

where {ni, ^2} are the chemical potentials and /liN^^'' is a sum over the two different gases. 
We then use a thermodynamic variational principle 

n{H)<n{H') + {{H-H%, (9) 



where 

(A), = 3 (10) 

and -ff* is a single particle trial Hamiltonian 

= E/ M^i'^I^.WlMa. + i^f I0.(r)p6l64- (H) 
This gives for the thermodynamic potential 

^{E) < n{El) + n{Ei) + (i^i - El)t + (^2 - H2)t + (^s)*. (12) 

The single particle states and (pa are to be determined such that they minimize the 
thermodynamic free energy Q. We now expand our field operators 

'4'i{r) =J2^»(^)"'» (13) 

a 

V'I(r)=E^:(r)4 (14) 

a 

Mr)=T.Mr)ba (15) 

a 

V^I(r)=EC(r)&l (16) 

a 

Inserting this expansion into Eq. |l] gives the Hamiltonians 



+ j drY, ¥'a(r)</'^(r)¥?^(r)(/?5(r)aj,aj3a^a5 

a/375 

E, = j rfrEC(r)[-^V2 + t/2(r)]0„(r)6t6^ 



(17) 



(18) 



E3 = w J drJ2 <(r)0K'^)'^7(i')'^'5('^)4^K"5- (19) 



We can now use the single particle Hamiltonians and rewrite the thermodynamic free energy 
as 

Q = + nl + J drn{r) (20) 

where 

= -i^/n[l-e-^(^«'-'^i)] (21) 

P a 

^2 = -4 E^^[l - e-^(^«'-'^^)] (22) 

P a 

(23) 

and Q = Qi + ^2 + ^3 with 



(24) 



+|Eb«(r)r(ArW(Ara) _ 1)), + ^, ^ ba(r)|>,(r)r(ArWArW), 



f^2(r) = EC(r)[-2;^V^ + C/2(r) - £;^)]0a(r)(7Vl^)), 



(25) 

+|El'^a(r)r(iVf (iVl^) - l)), + .2 E IMrmMrmNL'^NP), 



a 



Qsiv) ^wj:MrmMr)\'{Ni'^N^;\ (26) 

a/3 

Here we have used the independent particle properties of Hi and H2 with 

7V(^ = at,a„ (27) 

A^^^ = ^'l^'a (28) 

(alat,a,a,), = (iV«(iV«-l)), (29) 

(ttWt = (A^f (A^^^-l))t (30) 

(alG„6j,6^), = (iV«),(iv(^)), (31) 

(A^<i'^Ar(^^), = (A^<i'^)t(<^- (32) 



In order to calculate the single particle states and 0^ that minimizes the free energy fl 
we have to calculate the functional derivatives 

^ ■/rfrn«(r),0:(r)) = O (33) 
■JdrQ{^:{r),<p:{r)) = 0. (34) 
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For the Vli and parts we can use the results for the single condensate p!5|JT^ . From Eqs. 
(P3|) and ( piD we then get the equations 



2mi 



(35) 



+ EW^I^/5WP(A^^'^}^a(r) l'^/3WP(<^(A^^^^^a(r) = 

/3 P 



2m 



(36) 



+ EK'^)*l</>/3Wr(A^'^)a0.(r)+^El^/3(r)P(<)*(iV<'))*0a(r) = O. 



We also keep in mind the difference between the condensed phase and the normal phase 



concerning the single particle averages, see Ref. [13 



{N^^Nf), = (iV«),(iV«), a ^ ao,^ = 1,2 (37) 
{N^\n!^^ - 1)), = 2{N^^)l « ^ ao,^ = 1,2 (38) 

(A^S(^S - 1))* = No\n^^ - 1) ^ Nf. (39) 

We can now have three different regions depending on the parameters: Two condensates, 
a condensed phase and a normal phase or finally two normal phases. The three sets of 
equations are 

a) Two condensates 

The temperature is now below the critical temperature and T2 for both gases and we 
have 



where 





<^ao(r) 


=^i:Vao(r) 


(40) 








(41) 


1 r/iCr^l 1 V {'^[?{v\ 1 n^Ml 1 wT {v\ 
2m ^ ) ^ rn\^ ) ^ h'QK^ ) ) ^ ) 


rao 1^ ) 


— (r) 

— ^ao Veto \^ ) 




-^V^ + C/2(r) + ^2(2p^(r) + p^(r)) + wV^iv) 


(pair) 




(43) 








(44) 


r^(r-\ — A/"(l)l/o 






(45j 


pliy) = E {N^WU^)? 






(46) 


Po(r) = A^?|0ao(r)r 






(47) 






(48) 


r2(r) = <Vao(r)r+ E (^^'^l<^a(r)r 




(49) 


(7V(^)), = l/(e^(^"'-'^^) - 1) 






(50) 


(TVf), = i/(e/^(^^'^-''^) - 1). 






(51) 



b) A condensed phase and a normal phase 

The temperature is now below one of the critical temperatures. We give here the equa- 
tions with system (1) in the condensed phase, 



2mi 



2mi 



+ C/i(r) + vr{2pl{r) + p^(r)) + wr,{r 



2m' 



- V + U2{r) + 2v2p'{r) + wT2{r) 



<^ao(r) = ^i;Vao(r) 
99«(r) = Ei^Va(r) 
0«(r) = 0a(r) 



(52) 
(53) 
(54) 



with 



p'niy) = E (^P>*ka(r)p 



(55) 



p'{r) = Y.{NL%\Mr)\' 

a 
a 

r2(r) = iV«|y,,„(r)p + ^ (iVW),|^,(r) p. 
c) Two normal phases 

The temperature is here above the critical temperature for both gases. 



2mi 
2m-) 



+ Ui{y) + 2vip\r) + wTi{v) 
+ U2{v) + 2v2p\r)+wV2{r) 



0,(r) = i^f 0a(r) 



with 



a 

p\r) = j:{NL%\Mr)\' 

a 
a 
a 

The chemical potentials pi and 112 are calculated from the particle numbers 

a 

iV(2) = ^ 



^ g/3(i<;i^)-Mi) _ 1 

1 

r(2) 



(56) 
(57) 
(58) 
(59) 



(60) 
(61) 



(62) 
(63) 
(64) 
(65) 



In the condensed phase we have pi = E^'j and ^2 = ^an- 



(66) 
(67) 



III. THE NUMERICAL METHODS 



The numerical methods have been thouroghly explained in Ref. [|T3|. We are here only 
going to give the general ideas. The equations are highly nonlinear and therefore have to be 
solved iteratively. This means, we start by guessing the densities and solving the ordinary 
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differential equations. The resulting eigenvalues and eigenstates are then used to obtain 
new densities. This procedure is then iterated until we have a self consistent solution. In 
the experiments with Bose condensed gases the trapping potential has been approximately 
harmonic. The two gases have different external potentials because they are trapped in 
different spin states. This then gives us the two different external potentials 



Because of the harmonic external potentials, a natural approach would be to expand the 
solutions in the harmonic oscillator eigenfunctions 



The problem then reduces to an eigenvalue problem for the expansion coefficients a^. In 
the actual calculations we have at most four coupled equations that have to be solved 
simultaneously. This is not a problem since we have basically the same situation as in the 
single condensate case. The only drawback is the computing time which is doubled compared 
to the calculations in Ref. [|T^]. 

The method of expansion in harmonic oscillator eigenfunctions is especially good in 
asymmetric environments. However, if we have a spherically symmetric geometry, we can 
easily discretize our solutions on a lattice, ip{r) — > yjj, and the derivatives turn into differ- 
ences, which then gives us a tridiagonal eigenvalue problem whose solutions directly gives 
the desired quantities. This method works best in situations where the problem is effectively 
one dimensional. 



In this paper we calculate the density of the two-component condensate in a harmonic 
trap in two and three dimensions including asymmetry in the two dimensional case. We also 




(68) 



(69) 




(70) 



IV. RESULTS 
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present a phase diagram for a three dimensional spherically symmetric trap. The equations 
are put into dimensionless form by the scaling 

r = -r' (71) 

and X — Xx', y — Xy' with 

A = r- (72) 

[mim2] 2 [n^^n^^ny^ny^] s 



This gives us the dimensionless energies scaled to ^h{flxj^flx^fly^ ^2/2)^^^'\/^2/^i and 

A. The density calculations 

The most natural thing to calculate from the Hartree-Fock equations are the densities. 
In the spherically symmetric case we have the external potentials 

1 

1 



C/i(r) = -rrnny (73) 



U2{r) = ^m2fiy. (74) 

Throughout the calculations in this paper we have chosen the ratio between mi and m2 to 
be mi/m2 = 1 and the interactions vi and V2 to be vi = 0.02 and V2 = 0.01. In Fig.l we 
show the spherically symmetric density for the two condensates with N — 1000 particles in 
both gases and the interaction strength between the different particles is put to w — 0.012. 
Here we can see that the two different particles hardly know anything about each other. The 
inverse temperature is here /? = 0.3. In Fig.2 we have incresed the interaction strength w 
to w — 0.05. The condensate of the more weakly interacting particles is pushed away from 
the center of the trap and is forming a shell around the particles in the center of the trap. 
These results were obtained by spherically symmetric eigensolutions which means that we 
can not say anything about the existance of states with lower energy that could posess an 
asymmetric geometry. This suggests to look for a solution of an asymmetric two-component 
condensate. 
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It is extremely time and memory consuming to solve the full three dimensional asym- 
metric case. We therefore concentrate on the two dimensional two-component condensate. 
We then have the external potentials 

U,{x,y)^\rm{nlx' + niy) (75) 
U,(x,y)^^m2{nlx' + nly'). (76) 

In the two dimensional calculations we have chosen the asymmetry of the trap to be 

^k^VS i = l,2 (77) 

with Qai/^a2 — V^- III FigS- 3)-5) we show the density of the two condensates at the inverse 
temperature j3 = 1.0. If we start the iteration by putting the initial densities at for instance 
X — l,y — 1 and a; = l,y = — l,we end the iteration with two separated condensates aligned 
in the weaker trap direction. In Fig.3 we have used w — 0.1 and in Fig.4 we have incresed the 
interaction to w = 0.3, where we can see that the two condensates get further pushed away 
from each other with incresing interaction strength. However, if we start by putting both 
condensates at the center of the trap, we get a stable solution that is symmetric also in the 
weaker trap direction and docs not show a separation into two displaced condensates. This is 
shown in Fig. 5. The ground state energy for this case is higher than in the purely asymmetric 
situation which suggests that the spontaneous symmetry breaking occures. This is more 
clearly seen if we look at the two dimensional rotationally symmetric potential. Starting the 
iteration with the initial densities away from the center, gives us two condensates separated 
into two distinct peaks. This is shown in Fig.6 with /3 = 1.0 and w — 0.1. In Fig. 7 we have 
the same situation but now we start the iteration with the densities at the center of the 
trap. The calculation converges nicely and the stable solution shows two condensate peaks 
sitting on top of each other. In this case the interaction strength is not strong enough to 
create a ring of the lighter atoms around the center, which was seen in the three dimensional 
calculations in Fig. 2. 
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B. The phase diagram 



The Hartree-Fock equations do not describe exactly the transition region between an 
ordinary gas and a Bose condensed one, but they do suggest a general view of what is 
going on. If we define the critical temperature as the temperature where the number of 
particles in the condensate goes to zero, we can calculate this number of particles with the 
equations for (T < Tc) at different temperatures until we reach the point where all particles 
in the condensate have been depleted. In Fig. 8 we show three phase diagrams as function 
of temperature and interaction strength w, with a fixed number of particles A^*^^^ = 100 and 
varying between 100 and 728. We see that three regions need to be covered. We first 
solve the equations for two condensates with a fixed w and increse the temperature until one 
of the condensates disappears, where we find our phase transition. Incresing the temperature 
even more means that we are in the one-condensed-one-normal-phase region and we have 



to use Eqs. (p^)-(|5^). Incresing the temperature further, finally destroys the remaining 
condensate and we have found our second transition point. This kind of calculation has 
been repeted for different interaction strengths and particle numbers. Adjusting the particle 
numbers so that the two transition lines almost coincide, we find that there exists a region 
where one can go from two condensates to normal gases by lowering the temperature. This 
situation is shown in Fig. 8b. The phenomena exhibited in our phase diagrams may not 
give an accurat picture of the real systems, but they can be believed to suggest the trends 
expected in the actual experiments. 



V. CONCLUSIONS 

In this paper we have, for the first time, presented density calculations of a two- 
component condensate as functions of temperature. We have found a symmetry breaking 
in the two dimensional case, where the two-component condensate sitting at the center of 
the trap on top of each other posess a higher energy than the situation where the two 
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condensates have separated and form two individual peaks. This is a case of spontaneous 
symmetry breaking. We also present a phase diagram which describes the different regions 
with a two-component condensate, one single condensate and finally no condensate as func- 
tion of temperature and the interaction strength between the two different atoms. The phase 
diagram has been presented for three different pairs of particle numbers and calculated in 
three dimensions with spherical symmetric external potentials. This means that we can not 
see the symmetry breaking which we saw in the two dimensional calculations which, in fact, 
gave an asymmetric solution in the symmetrical situation. The Hartree-Fock equations do 
not describe the transition region well because they neglect all critical fluctuations. They 
do, however, give us an idea what may be happening around the transition point. At some 
critical particle number we find a region where it is possible by changing the temperature at 
fixed interaction strengths w, to go from a two-component condensate to two noncondensed 
gases. 

In these calculations, we have only used a few hundred particles. To increse the number 
of particles to realistic values (~ 10®) is beyond the capacity of available computers. So is 
a full treatment of totally asymmetric three dimensional traps. With this small number of 
particles, the interaction strength, w, needs to be chosen unrealistically large to bring out 
the observed features. The small number of particles may give rise to finite size effects as 
for instance regions with an inconsistant fraction of condensate particles calculated from the 
two-component condensate equations and the single condensate equations. With incresing 
particle number, we may expect the inconsistensy to disappear and give a crossing from 
a two-component condensate into a phase with two normal components for all interaction 
strengths w. 

The numerics was performed with a grid method in the three dimensional spherically 
symmetric case. This is a very stable and fast method. The only drawback is that it 
can in practice only be used in situations where the problem is effectively reduced to a one 
dimensional. The expansion method, on the other hand, is very fragile in the two dimensional 
calculations, and great care has to be taken in order not to get runaway iterations which do 
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not converge to physically stable solutions. 

Finally we want to speculate on some features observed in our numerical computations. 
They may be due only to shortcomings of the numerical approach or limitations of the HF 
method, but they point to interesting possibilities in the behaviour of the real systems. 

First we look at situations hke in Fig. 8b. There the two-condensate and no-condensate 
boundaries are very close to each other; they may, in fact, be found to cross. In such regions, 
we find that the critical temperatures found from the two-condensate and from the single- 
condensate sides are inconsistent. If we could trust the HF-calculations, this may indicate 
a hysteresis signalling a change to a first order transition for one of the components. No 
such conclusion can be proven from the HF-approach, but it points to the possibility that 
the change of one component can qualitatively affect the behaviour of the other one even to 
the extent that its order may change. 

Secondly, the use of too few states in the calculations gives a distorted and asymmetric 
solution. In two dimensions and for not too high temperatures (/3 — 1.0), we get smooth 
densities for about ten states in each direction. Because of the unrealistically large values for 
w used, the two condensates repell each other strongly. This shows up as a drift instability 
which eventually develops into an oscillational instability in the solution, if the computer 
iteration is continued well after a stable solution is found. It seems that this can be over- 
come by incresing the number of states involved. The physical contents of our numerical 
observation may be that it signals the break up of a solution only locally stable. Thus we 
find it to occur much more readily for the symmetric situation, whereas the asymmetric 
solutions are much more stable. 

We have thus found our numerics to indicate the change of order of the transition and 
the instability of some locally stable solutions owing to the interaction between the two 
condensates. Admittedly neither our numerical method nor our theoretical formulation 
(HF) allows any claims to the reahty of the effects. They do, however, offer challenging 
possibilities for further experimental and theoretical investigations. 

Unfortunately, it seems to be difficult to approach these problems from the two- 
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component Gross-Pitaevskii equations. As these omit the atoms above the condensate, 
it is difficult to compute the effects on one condensate by the properties of the other one. 
On the other hand, any phase transition theory superior to the HF-approach appears to 
offer unsurmountable computational difficulties. What progress can be achieved on these 
difficult questions remains to be seen. 
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FIGURE CAPTIONS 



Figure 1: The spherically symmetric density of the two condensates. The two condensates 
hardly see each other at w = 0.012. The insert shows the excited part. The temperature 
was here T — 3.33 and the interactions vi = 0.02 and V2 = 0.01 with A^*^^) = 1000 and 
A^(2) = 1000. 

Figure 2: The same situation as in Fig.l with a stronger mutual interaction w — 0.05. 
Here we can see that the A" -particles are pushed away from the center. 

Figure 3: The two dimensional asymmetric case were the condensates have formed two 
separated peaks around the center of the trap. The interaction is here w = 0.1 with Q^/^y = 
-\/8, 0,ai/^a2 — and the temperature T — 1.0. The iteration is here started away from 
the center with two Gaussian condensate densities at (1,1) and (1,-1)- The ground state 
energies are here ^^q^'* = 3.813 and E^^ = 2.924. 

Figure 4: The same situation as in Fig. 3 with the stronger interaction w — 0.3. The 
two condensates are pushed further away from each other because of their stronger repulsive 
interaction. 

Figure 5: The same situation as in Fig.3. The iteration is here started with the con- 
densates in the center of the trap with w — 0.1 and T — 1.0. The condensates sit on top 
of each other. The ground state energy is higher than in the asymmetric formation, with 
E^^^ =4.114 and E^^^ = 3.150; cf the values given in Fig. 3. 

Figure 6: The two dimensional spherically symmetric situation were we find two sepa- 
rated density peaks with the ground state energies E^^^ — 3.80, E^^^ — 2.97 and w — 0.1. 
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The different lines sfiow tfie density contours as indicated on tlie left side of the figure. The 
iteration is started away from the center. 

Figure 7: The two condensates situated in the center of the trap on top of each other. 
The iteration is started with two Gaussian densities positioned at the center of the trap. The 
ground state energy in this case is found to be slightly higher than the off-center iteration 
in Fig.6. The energies are here E^^^ = 3.81 and E^^^ = 2.98. 

Figure 8: The phase diagram described for three different pairs of particle number. In 
a) we have A^(^) = 100 and Ar(2) = 100. Here we can see that the transition line between 
region II and III (one condensate and no condensate) is not sensitive to changes in w. In 
region II the A?"^^^ -particles are condensed and the A^^^^ -particles are in the normal phase. 
In b) we have = 100 and A^^^^ = 617, and we find a region where we have a transition 
between a two-component condensate (I) and no condensate (III) . In c) we have incresed the 
relative particle number to N^^^ — 100 and N^'^^ — 728. Region I with the two-component 
condensate is strongly supressed as a function of w, whereas the transition between one 
condensate and no condensate (II and III) is not sensitive to changes in w. In this case, for 
region II, the A'"^^^ -particles are condensed and the A'"^^^ -particles are in the normal phase 
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Figure 3. P. Ohberg and S. Stenholm 
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Figure 7. P. Ohberg and S. Stenholm 
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